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Abstract. Wc consider the stationary state properties of the reduced density matrix as well 
as spin-spin correlation functions after a sudden quantum quench of the magnetic field in the 
transverse field Ising chain. Wc demonstrate that stationary state properties are described by a 
generalized Gibbs ensemble. We discuss the approach to the stationary state at late times. 
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1. Introduction 

This is the second of two papers on the quench dynamics of the transverse field Ising chain. In the 
first part of our work [1], which in the following we will refer to as "paper I" , we focussed on the time 
dependence of the longitudinal spin correlations. The present manuscript gives a detailed account 
of properties in the stationary state at infinite times after the quench. An important motivation 
for studying problems of nonequilibrium time evolution in isolated quantum systems is provided by 
recent experiments on trapped ultra-cold atomic gases [2, 3, 4, 5, 6, 7]. These experiments suggest 
that observables such as multi-point correlation functions generically relax to time independent 
values. Such a behaviour at first appears quite surprising, because unitary time evolution maintains 
the system in a pure state at all times. The resolution of this apparent paradox is that in the 
thermodynamic limit, (finite) subsystems can and do display correlations characteristic of a mixed 
state, namely the one obtained by tracing out the degrees of freedom outside the subsystem itself. 
In physical terms this means that the system acts as its own bath. An important question is how to 
characterize the reduced density matrix describing the stationary behaviour. Intuitively one may 
expect that conservation laws will play an important role, which in turn poses the question whether 
quantum integrable systems exhibit qualitatively different stationary behaviours when compared 
to generic, nonintegrable ones. These issues have been addressed in a number of recent works 
[8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 
34, 35, 36, 37, 38, 39, 40, 41, 42]. Most of these studies are compatible with the widely held belief 
(see e.g. [8] for a comprehensive summary) that the reduced density matrix of any finite subsystem 
(which determines correlation functions of all local observables within the subsystem) of an infinite 
system can be described in terms of either an effective thermal (Gibbs) distribution or a so-called 
generalized Gibbs ensemble (GGE) [9]. It has been conjectured that the latter arises for integrable 
models, while the former is obtained for generic systems. Evidence supporting this view has been 
obtained in a number of examples [9, 10, 12, 13, 14, 15, 16, 19, 21, 22, 23]. On the other hand, 
several numerical studies [24, 26, 27, 37, 38] suggest that the full picture may well be more complex. 
Moreover, open questions remain even with regard to the very existence of stationary states. For 
example, the order parameter of certain mean-field models have been shown to display persistent 
oscillations [43, 44, 45, 46, 47, 48]. Non-decaying oscillations have also been observed numerically 
[26] in some non-integrable one-dimensional systems. This has given rise to the concept of "weak 
thermalization" , which refers to a situation where only time-averaged quantities are thermal. 

We will show in this manuscript that for a quench of the magnetic field in the transverse 
field Ising chain, the reduced density matrix is described by a GGE. This establishes that any 
local observable is described by the GGE. However we will also show that while some two-point 
observables approach their asymptotic values relatively quickly, others will do so only after times 
exponentially large in the separation between the two points. In practice this precludes experimental 
or numerical detection of stationary behaviour for these observables. 

1.1. Quench protocol and observables 

In the following we focus on a global quantum quench of the magnetic field in the Ising Hamiltonian 

L 

H(h) = -jY,[<r*<r* + i+k*]], (1) 

3 = 1 

where af are the Pauli matrices at site j, we assume J > 0, h > and we impose periodic 
boundary conditions = a". We assume that the many-body system is prepared in the ground 
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state l^o) °f Hamiltonian H(ho). At time t = the field h is changed instantaneously to a different 
value h and one then considers the unitary time evolution of the system characterized by the new 
Hamiltonian H(h), i.e. the initial state l^o) evolves as 

|*o(*))=e- itff W|* ). (2) 
The above protocol corresponds to an experimental situation [2, 4, 5], in which a system parameter 
has been changed on a time scale that is small compared to all characteristic time scales present in 
the system. 

In equilibrium at zero temperature (and in the thermodynamic limit) the Ising model (I) 
exhibits ferromagnetic (h < 1) and paramagnetic (h > I) phases, separated by a quantum critical 
point at h c = 1. The order parameter for the corresponding quantum phase transition is the 
ground state expectation value (<rj). The Hamiltonian (1) can be diagonalized by a Jordan- 
Wigncr transformation, which maps the model to spinless fermions with local annihilation operators 
Cj, followed by a Fourier transform and, finally, a Bogoliubov transformation. In terms of the 
momentum space Bogoliubov fermions a k the Hamiltonian is diagonal: 

H{h) = ^e h {k)a\a k , e h {k) = 2J^l + h? - 2/icos(fe). (3) 

fc 

Details and precise definitions are given in Appendix A of paper I. 

In the following we focus on the one and two-point functions of the order parameter 

P [t) - <*o(t)|*o(*)> ' [ ) 

pXX ^ t] = ^mmm)^ ' pTm = ^ - {pX[t)f ' (5) 

the transverse spin correlators 

z(t) = (MQWW)) zz(e t) = (Mt)\*3 + tf\Mt)) 2 

9 W (*o(*)|*o(*)> ' Pc { , ) (*o(*)|*o(t)> {P ( )) ' [ ) 

and the reduced density matrix of subsystem A 

p A (t) = Tt A |* (*)X*o(*)|. (7) 
Here A U A is the entire system and in practice we will take A to consist of I consecutive sites along 
the chain. Tr^ denotes a trace in the space of states describing the lattice sites in A. 

Following paper I, we divide the time evolution of a two-point function for a fixed distance 
t between the operator insertions into three regimes, which are determined by the propagation 
velocity v(k) — rfg ^ fc - ) of elementary excitations of the post-quench Hamiltonian. For a given final 
magnetic field h, the maximal propagation velocity is 

Wmax = max \e' h (k)\ = 2 J min[/i, 1] . (8) 

The three different regimes are: 

• Short-times v max i <C £■ 

• Intermediate times v max t <~ I. This regime is of particular importance for both experiments and 
numerical computations. A convenient way of describing this regime is to consider evolution 
along a particular "ray" k£ = w max i in space-time, see Fig. 1. In order to obtain an accurate 
description of the dynamics at a particular point along this ray, one may then construct an 
asymptotic expansion in the single variable I around the space-time scaling limit v max t, i — > oo, 
k fixed. 




Figure 1. Left panel: for intermediate times Vm^t ~ £ the behaviour of p aa (£,t) is most 
conveniently determined by considering its asymptotic expansion around infinity ("space-time 
scaling limit") along the ray u m axt = k£. This viewpoint is appropriate for any large, finite t 
and i. Right panel: the asymptotic late-time regime is reached by considering time evolution at 
fixed t. To describe this regime one should consider an asymptotic expansion of p aa (£, i) around 
t = oo at fixed £. 



• Late times v max t I. This includes the limit t — > oo at fixed but large £. In this regime it is no 
longer convenient to consider evolution along a particular ray in space-time. In order to obtain 
accurate results for the late time dynamics, one should construct an asymptotic expansion in 
t around infinity, see Fig. 1. 

It is important to note that in general taking k — > oo in the space-time scaling limit docs not 
necessarily reproduce the late time behaviour at fixed, asymptotically large t. In other words, in 
general we have 

lim' W p° a (£,t)^Bm' lim' p™(l,t), (9) 

K fixed 

where lim' denotes the leading term in an asymptotic expansion around the limiting point. We 
will show that the limits do commute for p xx (£,t) for most but not all quenches, but they never 
commute for p z c z (£, t) . 



1.2. Longitudinal versus Transverse Correlators 

A global quantum quench of the transverse field in the Ising model is an essentially ideal testing 
ground for many ideas related to thermalization [49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61]. 
An example is the issue to what extent the late time behaviour after a quench may depend on the 
observable under consideration. In this regard, one may expect the locality of observables relative to 
the elementary excitations of the model to play a role [49, 26, 50]. In the Ising chain, the transverse 
spin operator ct| is local with respect to the fermionic degrees of freedom Cj, while the order 
parameter erj is in general non-local. The non-locality makes the longitudinal correlators difficult 
to analyze and general analytic results have been reported only recently in our short communication 
[51] (the stationary properties in the special cases ho — 0, oo were obtained before in Ref. [55]). The 
results reported in [51] for the correlation length in the stationary state were found to be described 
by an appropriately defined GGE. 

In contrast to order parameter correlators, the time evolution of one and two point functions 
of the transverse spins is straightforward to analyze (as the latter are fermion bilinears) and 
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has been known since 1970 [52]. In spite of its simplicity, the connected two-point function 
p zz (£, t) (corresponding to density correlations in gases) exhibits an interesting and quite general 
phenomenon, namely a cross-over in the relaxational behaviour at an exponentially large time scale 
Across ~ e £ - As a consequence the truly stationary regime of p zz (£,t) is observed only at very late 
times t > Across, which in general is a serious limitation. 



1.3. The quench variables 

As shown in Appendix A of paper I, both the initial and final Hamiltonians can be diagonalizcd 
by combined Jordan- Wigncr and Bogoliubov transformations with Bogoliubov angles 9® and 9k 
respectively 

e»* = , H - etk , e< = - . (10) 

V 1 + h 2 - 2h cos k yjl + hi - 2h cos k 

The corresponding Bogoliubov fcrmions are related by a linear transformation characterized by the 
difference A& = 9k — in order to parametrize the quench it is useful to introduce the quantity 

hh - (h + h ) cosfc + 1 
cosA fc = — , — ; (11 
yjl + h 2 - 2/icos(fc) A /l + h%- 2h cos(k) 

We note that cos A^ is invariant under the two transformations 

(h ,h) ^ (h,h ) and {h ,h)^ (^l) ■ (12) 
However, we stress that the quantum quench itself is not invariant under the maps (12). 

1.4. Generalized Gibbs ensemble 

The density matrix of a GGE can be cast in the form [62] 

Pgge= -J— e-^^'™, (13) 

^GGE 

where I m are some integrals of motion and Zqge ensures the normalization condition TrpcGE = 1- 
In Ref. [9], Rigol, Dunjko, Yurovsky and Olshanii proposed that an integrable system after a 
quantum quench in the infinite time limit is in fact described by a GGE, where the I m represent 
a complete set of independent integrals of motion and the Lagrange multipliers A m are fully 
determined by the initial state |^o) through the conditions 

Tr[p GGE / m ] = (*o|/m|*o> • (14) 

As we are particularly interested in the thermodynamic limit, we will employ the following, 
somewhat narrower, definition of a GGE for the description of stationary properties after quantum 
quenches. 

• First, as the entire system will be in a pure state at all times, it cannot be described by a 
density matrix corresponding to a mixed state. We therefore take Pgge to be the reduced 
density matrix of an infinitely large subsystem A in the thermodynamic limit. From this 
reduced density matrix all multipoint correlation functions can be obtained. Moreover, one 
can extract the reduced density matrix of any finite subsystem A\ by [15, 16] 

p Al (t = 00) = Tr Ai p GGE . (15) 
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Considering an infinitely large subsystem A has the advantage that the density matrix 
defining the GGE can be expressed in a simple way even for more complicated, interacting 
intcgrablc models, whereas the corresponding formulation for a finite subsystem becomes 
difficult. However, in cases where the subsystem A x consists of several disjoint blocks, the 
reduced density matrix will not have a simple form like (16) even in the absence of interactions 
[63]. 

The precise sequence of limits we have in mind is the following. The entire system is decomposed 
in a subsystem A and its complement A. We then take the thermodynamic limit, keeping A 
fixed. Finally we take the limit of A becoming infinite itself. The operator pgge is the reduced 
density matrix of subsystem A in this limit. 

The importance of considering the reduced density matrix when applying GGE ideas has been 
emphasized previously in Rcfs. [15, 16, 19, 26, 27]. 

• Another crucial point, that perhaps has not yet received the attention it deserves, concerns the 
issue of what integrals of motion I m ought to be included in the definition (14) of the GGE 
density matrix. Indeed, any quantum system, integrable or not, has many integrals of motion. 
For example the projectors I n = \ip n ) (tp n \ on Hamiltonian cigenstates H\i[) n } — E n \i[) n ) are 
integrals of motion [H, I n ] = [I m , I n ] = 0. Clearly, such conservation laws cannot play any role 
in determining the late time behaviour after a quantum quench. Following [64] we therefore 
propose to include only local integrals of motion in (14). These are characterized by arising 
from an integral (or a sum in the case of lattice models) of a local current density J n {x) as 
I n = J dxJ n (x), in the same spirit of the well known Noether theorem in quantum field theory. 
This locality requirement is similar in spirit to the additivity requirement of Ref . [65] . 

With these remarks in mind, the generalized Gibbs ensemble for the Ising chain is then defined 
by the density matrix [9] 



1 / x 

PGGE = 7} ex P 

^GGE 



Pk£ka\a k , (16) 



k J 



and the Lagrange multiplier f3 k are fixed by the initial state |*o) through the equations 

(*ol4 a fcl*o> = Tr[p GGE a[a fe ] . (17) 

For a quench of the magnetic field in the Ising chain, these are easily solved with the results 

fik£k = 2arctanh(cos Afe) . (18) 

However, the integrals of motion n k = ot\ctk use d in (16) arc non-local in space. A priori this is a 
serious problem as the integrals of motion defining a GGE must be local, as we pointed out above. 
However, as shown in Ref. [64], the prescription (16) works, because in the particular case of the 
Ising model the local integrals of motion can be expressed as complicated linear combinations of 
nk = u\ctk- Hence for an infinite subsystem (16) is equivalent to a GGE defined using the local 
integrals of motion. 



1.5. Organization of the manuscript. 

The present manuscript is organized as follows. In section 2 we give a detailed summary of our 
main results. In Sec. 3 we show that the reduced density matrix is described by a GGE. In Sec. 4 
we compute the stationary values of the longitudinal two-point correlation function, while in Sec. 
5 we study transverse correlations and their approach to the GGE. In Appendix A we review some 
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mathematical theorems we need to evaluate the asymptotic behaviour of longitudinal correlation 
functions. 

2. Summary of Results 

2.1. Reduced Density Matrix in the Stationary State 

Given that we start the system out in a pure state |^o) , its full density matrix is the one-dimensional 
projector 

p(t) = |* (t)X*o(t)|. (19) 

The reduced density matrix of a subsystem A is then obtained as 

p A (t) = Tv A (p(t)), (20) 

where A is the complement of A. We have obtained the following result for p A (t = oo): in the 
thermodynamic limit, for the case where A is itself infinite, the reduced density matrix in the 
stationary state is equivalent to that of the generalized Gibbs ensemble (16) 

Pa{1 = oo) = pgge- (21) 

In particular, this implies that arbitrary local multi-point spin correlation functions within 
subsystem A can be evaluated as averages within the GGE. However, the calculation of such averages 
represents itself a formidable problem. We have determined these averages for the most important 
cases of the two-point functions p xx (£,t = oo) and p zz (£,t = oo) and present explicit expressions 
below. We stress that although a x is not local in terms of fermions, correlation functions involving 
a x can nevertheless be calculated from pgge- For correlation functions involving an even number 
of a x, s this is because the Jordan- Wigner strings cancel and one is left with an expression that 
involves only fermions within subsystem A. Correlation functions involving an odd number of a x, s 
vanish as a result of the Z 2 symmetry, which is restored in the GGE. 

2.2. Longitudinal Correlators in the Stationary State 

We find that the infinite time limit of the two-point function lim^oo p xx (£, t) for fixed but 
asymptotically large £ is given by 

p xx (£~> l,t = oo) = C x (£)e~ e ^ [l + o(£ )] . (22) 

Here £ and C x (£) are functions of the initial (ho) and final (h) magnetic fields. The inverse 
correlation length is 

r 1 = H (h - l)e H (h - 1) In [min(/i , hi)] - ln[x+ + x_ + 8 H ((h - l)(h - l))y/4x+X-] , (23) 
where 9h{x) is the Heaviside step function and 

[min^fo-^inmintVV 1 ^!] 1 + hh + y/(h? - l)(h 2 - 1) 

X± = 4 ' kl = hTho ' (24) 

The function C x (£) describes the subleading large-f asymptotics and takes the following form 
(i) Quench within the ferromagnetic phase (ho, h < 1). 

C x (£) = l^WOE^l) s c x FF . (25) 
V ' 2^1-hho^/T^hJ FF V ' 
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(ii) Quench from the ferromagnetic to the paramagnetic phase (ho < 1 < h). 



C x {t) 



'hy/T 



h 2 
a 



C FP . 



h + h 

(iii) Quench from the paramagnetic to the ferromagnetic phase (ho > 1 > h). 



C x (£) 



h0 k cos(^arctan V^IEgMES ) s CpF W . 



1 + h h 



(iv) Quench within the paramagnetic phase (1 < ho,h). 
C x (£) = C£p(£) = < 



h Vh(hh -l + y?(h 2 -l)(hl-l)) 

4V¥(feg-l)3/4(fe fe-l)3/2 (;i _ feo) 

h(hg-h)y/hl-\ 
(h+hg)(hhg-l) 



£- 3 / 2 if 1< h < h , 
if 1 < ft< fto- 



(26) 



(27) 



(28) 



Here arctan(|x|) € [0,7r/2). We see that in most cases C x (f) tends to a constant value at large I. 
The exceptions are quenches from the paramagnetic to the ferromagnetic phase, where C x (£) tends 
to an oscillatory function with constant amplitude, and quenches to a larger magnetic field within 
the paramagnetic phase, where C x (t) decays like a power law with exponent —3/2. For the special 
cases ho = and h = oo the asymptotic behaviour of p xx (£ 1, t = oo) agrees with the results of 
Rcf. [55], which were obtained by different methods. 

2.3. Transverse Correlators in the Stationary State 

The stationary behaviour of the connected longitudinal two-point function is 

l,t = oo) ~ CT'e-^lHOfr 1 )), (29) 
where the transverse correlation length £ z , the exponent a z and the amplitude C z are given by 
C 1 = |Inft | +min(|ln/i |,|ln/i|), 

1 if | Inft| > | lnfto| , 

n = < if h = l/h , 

1/2 if |Inft| < |lnft | , 



(30) 



C z = < 



\h -l/h \ hg — h 
47T hhg — 1 
(h-l/h) 2 

27T 

(h-l/h)y/\ho-l/hg\(h -h) 



hg-h 



sgn(lnfe)|ln hg\/2 



hg{hhg — l) . 



lnh| + |lnh | 



if |ln/i| > |ln/i | , 
if ho ~ l/h , 

if |lnft| < |ln/i |. 



(31) 



2.^. 77ow Zon^ one has to wait before one can "see" the GGE 

Given that the infinite time behaviour of the reduced density matrix is described by a GGE, a 
natural question to ask is how long we need to wait in order to be able to detect the convergence 
of a given observable to its stationary value. Clearly, to answer this question we require knowledge 
of the full time evolution of correlation functions and not only their stationary values. This is very 
difficult in general, but can be studied in detail for the simple yet instructive case of the transverse 
spin-spin correlation function. 
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2.4-1- Transverse Correlations The connected transverse spin-spin correlation function has the 
following integral representation 

2 



- Y[ e 9kj cos A fe . - i sin A fcj cos (2e h (kj)t) \ - (32) 
j=i J 
The integral representation can now be used to obtain asymptotic expansions for large £ and t in 
the two limits of interest: 

(i) In the space-time scaling limit £, ti max t — > oo with v max t/£ fixed, the asymptotic behaviour 
can be evaluated by means of a stationary phase approximation. This shows that the leading 
behaviour is aT 1 power-law decay, while sublcading corrections are power laws as well, i.e. 

p?(t= V ^,t)~ D lW + o(t-i), (33) 

where D z (t) is sum of a constant contribution and oscillatory terms with constant amplitudes. 

(ii) In the late time regime at fixed, large £, we find that p zz (I, t) decays as a power law in t to a 
stationary value that is exponentially small in £ 

p zz (£,t) ~ pf (*,oo) + E ' {t) £~ m ' +o{t-*'*) , (34) 

where E z {t) = ^2 q=Q K A q cos (2teh(q) + <p q ) with constants A q , ip q , and the stationary value 
p zz (£, oo) is given above in (29). Crucially, p zz (£, oo) oc e _£ /^ is exponentially small in £. The 
inverse correlation length £j x is given by 

C 1 = min(| log/iol, I log/i|) < C 1 ■ (35) 

We note that in the space-time scaling limit exponentially small terms such as p zz (£, oo) will always 
be negligible compared to the dominant power law behaviour. On the other hand, in the late time 
regime there exists a cross-over time scale 

after which the stationary behaviour becomes apparent. Importantly, this time scale is exponentially 
large in the separation £. As we have alluded to before, the space-time scaling limit is a convenient 
way of obtaining the behaviour of p zz (£,t) for general large £ and t. Contributions that are 
exponentially small in £ arc not included in the asymptotic expansion obtained in the space-time 
scaling limit. This observation allows us to define a cross-over time scale t\ between the intermediate 
and late-time regimes by the requirement that the leading asymptotics in the space-time scaling limit 
(oc £ 2 /t 3 ) becomes comparable to the neglected exponentially small terms 0(e~ e ^ z , e~ e ^ z /t 3 / 2 ). 
This gives 

Jt\~e 21 ^. (37) 

This estimate suggests that the result obtained in the space-time scaling limit will give a good 
description of p zz (£, t) for fixed separation £ up to times of order t\. When | log ho\ < \ log h\, the 
two time scales are comparable t* ~ t* 2 . This means that in practice the space-time scaling limit 
provides a good approximation for p zz (£, t) all the way up to the stationary regime (see Fig. 2). On 



Quantum Quench in the Transverse Field Ising chain II 



10 




Figure 2. Absolute value of the connected transverse correlator p zz (£,t) — p zz (£,t = oo) as 
function of time at fixed £ = 10. Initial and final magnetic fields ho = 0.6 — > h = 0.5 are chosen 
such that a crossover between a t~ 3 and the asymptotic t~ 3 / 2 regime occurs at accessible times. 



the other hand, when | log/io| > | log/i| there is a large time window t* 2 > t > t\ in which a t~ 3 / 2 
power-law decay can be observed. 

The implications of these considerations for experiments or numerical calculations on finite-size 
systems are as follows: 

• First and foremost, in order to observe the GGE either numerically or experimentally, it is 
essential to focus on a suitable observable. For the case of the TFIC this would be e.g. the 
one-point function p z (t). 

• In a finite system of size L, boundary effects (such as reflection from the boundaries or 
completing a full system traverse in the periodic case) will become important at a timescale 

tf a ~L/v maK . (38) 

In order to observe stationary behaviour, i.e. the GGE, this timescale should be larger than 
the cross-over scale t* 2 

tf s > t*. (39) 

In practice this requires extremely large system sizes. Just to quote some numbers: if we could 
simulate chains of length L = 10000 (which is at least one order of magnitude larger than 
what can be achieved with currently available numerical algorithms), the GGE can be "seen" 
at best for correlations evaluated at distances of about 10 lattice sites. Correlation functions 
evaluated at larger distances may appear stationary, but this behaviour will not be related to 
the GGE. 

• The behaviour for large v max t < u max ifs and distances I (compared to the lattice spacing) for 
a finite system will generically be described by the intermediate time regime and in particular 
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the result obtained in the space-time scaling limit. Time scales necessary to observe stationary 
behaviour will generically be inaccessible. 



2.4-2. Longitudinal Correlations As a result of their non-local nature with respect to the Jordan- 
Wigner fermions these are much more difficult to determine than the transverse correlations. 
Detailed results for the time evolution of one and two-point functions of are presented in paper 
I. The implications of these results for the present discussion are summarized as follows. 

• For quenches originating in the ferromagnetic phase, the stationary value of p xx (£,t) emerges 
after a time ijj that scales as a power-law with the spatial separation v max t^ ~ £ 4 ^ 3 , cf. Eq. 
(28) of paper I. Hence it is straightforward to see the approach (as a function of time) of the 
two-point function to its stationary values, which are given by the GGE. 

• For quenches within the paramagnetic phase p xx (£,t) exhibits an oscillatory power- law decay 
in time towards its stationary value, which is exponentially small in I. Hence, in complete 
analogy to the case of the transverse two-point function, the time scale t\ after which the 
stationary behaviour reveals itself, is exponentially large 

t* 2 oc e 2 ^ 3 «, (40) 

and very difficult to observe in practice. 



2.5. Generalized Gibbs Ensemble and Conformal Field Theory. 

In Refs. [11, 12] it has been shown that for conformal field theories the long time limit of correlation 
functions after a quantum quench described by a conformally invariant boundary state is described 
by a thermal (Gibbs) ensemble. Since CFTs are prototypical integrable field theories, this result 
apparently contradicts the general expectation that integrable models should be described by 
appropriately defined GGEs. We now show that for the case of the Ising field theory there is 
in fact no contradiction and that the peculiarity of the CFT result can be traced back to the 
particular (non-generic) initial state considered in Refs. [11, 12]. 

As discussed in paper I, CFT describes the scaling limit of the Ising model when the final gap 
is send to zero. The scaling limit of the transverse field Ising chain is (a is the lattice spacing) 

J^oo, /i -> 1 , a ^ 0, (41) 
while keeping fixed both the gap A and the velocity v 

2J|l-/i|=A, 2Ja = v. (42) 
In this limit the dispersion and Bogoliubov angle become 

s(q) = y/A 2 + v 2 q 2 , 9 h (q) -> sgn(ft - l)arctan(^) . (43) 

Here the physical momentum is defined as q = fc/a , with — oo < q < oo. In our quench problem 
both the initial and the final magnetic field are scaled to the critical point, i.e. we need to take 

ft ->l, 2J|1 - /i 1 = A = fixed. (44) 
Thus in the scaling limit we have 

cosA fc ^ - AA + (vg) 2 

^(A 2 + (vq) 2 )(A 2 + (vq) 2 ) 
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For quenches to the quantum critical point, which is described by the Ising conformal field theory 
with central charge c = 1/2, we need to furthermore take A — > 0, which results in 

cos A fe -> - " |g| . (46) 

In this limit the expression (18) for the Lagrange multipliers j3k becomes 

Pav\q\ = 2 arctanh— ^"^^ . (47) 

This shows that one obtains mode-dependent temperatures and hence a non-trivial GGE even in 
the case of a quench to the Ising CFT. 

However, the particular boundary states considered as initial states in Refs. [11, 12] correspond 
to the limit of infinitesimal correlation length, i.e. a very large initial gap Ao — > oo. In the limit 
v\q\/A -> 0, (47) becomes 

2 

Ao ' 

i.e. all mode dependent temperatures become equal, resulting in a thermal state. Thus the findings 
[11, 12] should not be interpreted as showing that CFTs thermalizc after quantum quenches. We 
have been informed by John Cardy that by perturbing the conformal boundary condition it is 
possible to show the emergence of a mode dependent temperature directly in CFT. 



P q = ~, (48) 



3. The infinite time limit and the generalized Gibbs ensemble 



In this section we consider the infinite time limit of the reduced density matrix pa of a subsystem 
A composed of I contiguous spins. To analyze pa, we first consider its building blocks, i.e. the 
two-point real-space correlation functions of fermions. It is convenient to replace the Jordan- Wigner 
fermions cj (as defined in Appendix A of paper I) by the (real) Majorana fermions 



a* = c]+ Cj 



i y — a( 



= i(c}- Cj ) 



which satisfy the algebra {of, a%} = 2Si n , {of, a^} = 2<5j„, {af , a^} 
fermions, the operator trj has the nonlocal representation 



while af is local 



of = ia^af . 
The correlation matrix T is defined as 

r r_j ••• rv 



(49) 

0. In terms of these Majorana 

(50) 

(51) 



r = 



r 



-ft 

-9-i 



91 
fi 



(52) 
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where the matrix elements are the fermionic correlators 
9n = ~ (afaf +n ) , 

f n + S n0 = (afaf +n ) = (a y a y ) Ml. (53) 

The matrix T is a block Toeplitz matrix, because its constituent 2x2 blocks depend only on the 
difference between row and column indices. It is customary to introduce the (block) symbol of the 
matrix T as follows 



r(k) = 



-/(*) 9(k) 

-g(-k) f(k) 



oo 

E 



— ink-p 

1 7) 



The functions f(k) and g(k) are 

f(k) = sin A k sm(2e h (k)t), 



g(k) = -ie 



(54) 



(55) 



cos Afe — isin Afe cos(2£(j(fc)i) , 

where e l6k and cos A/; are defined in (10) and (11) respectively, while 
. (h- h ) sink 

sm A fe = — . (56) 

VI + h 2 — 2h cos k\J 1 + Kq — 2h cosfc 

The infinite time limit of the fcrmion correlators (53) is straightforwardly obtained by Fourier 
transforming the functions f(k) and g(k) 



/f = (a?a? +J ) 



Sjo = 0, 



t— oo 



/7T 
-7T 



dfc 
27 



e -i0k e -i6 k cosAfc . 



(57) 



We can now use the Wick theorem to construct all correlation functions in the Ising chain!. 
Moreover, as shown in Refs [66, 67], the matrix Y determines the full reduced density matrix 
of the block A of t contiguous fermions (and hence spins, because contiguous spins are mapped to 
contiguous fermions, see Eq. (50)) in the chain 



PA 



1 21 

h E (IN 



2« 



where a 2n = <, a 2n -i 



w =0,l 1=1 

al, and 




W 

tanh — = T . 



(58) 



(59) 



Given pa one can calculate any local correlation function with support in A. 

In order to understand the properties of the stationary reduced density matrix, we should study 
the structure of the fermion correlations in Eq. (57). These resemble the analogous correlations of 
an Ising chain in equilibrium at finite temperature calculated in Ref. [52] 

fG8) 



dk 

2^ 



e -ijk e -i0 k tanh 



(60) 



J For t — ¥ oo, the expectation values of odd (with respect to fermion parity) operators vanish for any quench. Thus 
all nonzero correlation functions can be written in terms of two-point fermion correlators only. 
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where ((C)) denotes the equilibrium expectation value of the operator O at temperature 1//3. 
Comparing (60) to (57) suggests that if it were possible to find a (3 such that cos = tanh(/3efc/2) 
for any k, local properties of the (sub)systcm could be described by a thermal state at inverse 
temperature /3 _1 . It is easy to see that this is not possible: on the one hand, the dispersion relation 
is a continuous bounded function such that tanh(/3efc/2) < 1 for any finite f3. On the other hand 
we have | cos A | = | cos A„-| = 1, which implies that the equality cannot be fulfilled at least in some 
neighbourhood of the momenta and ir. 

Having ruled out the possibility of a thermal reduced density matrix in the stationary state, we 
next investigate the conjecture put forward by Rigol et. al. [9], according to which the stationary 
state of integrable models are described by a generalized Gibbs ensemble. The latter is obtained by 
maximizing the entropy, while keeping the energy as well as all higher conservation laws fixed. The 
resulting stationary reduced density matrix is of the form 

PGGE = T^—e-^ x ^ 1 (61) 

^GGE 

where the I m represent a complete set of independent, local integrals of motion and the Lagrange 
multipliers A m arc fully determined by the initial state \^o) through the conditions 

Tr[p GGE / m ] = (*o|/m|*o) • (62) 

We stress that locality of I m is an essential feature of this formulation. For the case of the 
Ising model, it is shown in Ref. [64] that the local integrals of motion can be expressed as linear 
combinations of the Bogoliubov fermion number operators = a\a^. Hence for the Ising model 
(and by analogy other theories with free fermionic spectra) one may replace the local charges I m in 
(61) by the number operators n/., despite the fact that the latter are non-local as discussed in the 
introduction. Thus, taking the Lagrange multipliers as A fe = /? fe e fe in order to facilitate the eventual 
identification of /3k as a mode-dependent temperature, the GGE takes the form 

Pgge= -^—e-^^i"". (63) 

^GGE 

The Lagrange multipliers 0k are fixed by the conditions 

Tr[p GGE n k ] = (*o|n fc |tfo) , (64) 
which can be solved in closed form by 

Pk£k = 2arctanh(cos A fe ) . (65) 

The fermionic two-point correlators in the GGE can be straightforwardly evaluated following the 
finite-temperature calculation and are given by 

gf GE) = -i f 2 e ~ iike ~ i9k tanh (^?) ' (66) 

The crucial point now is that the fermionic two-point functions g°° and g^ GE ^ are exactly the 
same, given that the Lagrange multipliers /3k ar e related to A& by (65), i.e. 

9T=9f GE) - (67) 

As the fermion correlators f°° and g°° completely fix the reduced density matrix pa (the matrix 
Wim in (58) can be determined by simply evaluating all fermionic two-point functions using Wick's 
theorem) , the equality (67) is lifted to the level of the reduced density matrices 

PA = Tr A (p GGE ). (68) 
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A characteristic feature of the GGE is that it retains detailed information about the initial state. 
This is in contrast to a thermal ensemble, which would correspond to fik = ft for any k and would 
depend on the initial state only through the average value of the energy. Given the structure of 
A fe there is no quench for which the f3 k is independent on k and hence the stationary state is never 
thermal. However, it is possible for the differences between GGE and thermal ensemble to be small 
for certain observables [49]. 

To summarize the results of this section, we have shown that the stationary reduced density 
matrix after an arbitrary quench of the bulk magnetic field in the Ising chain is equivalent to a 
generalized Gibbs ensemble. Hence all local correlation functions can be deduced from the GGE. 
The calculation of observables in the framework of the GGE generically remains a difficult problem. 
In the following sections we determine the stationary behaviour of some of the most important 
observables, namely the longitudinal and transverse correlation functions. 



4. Stationary value of the longitudinal correlation function 

The two-point function of a x is the expectation value of a string of Majorana fermions 

p xx {£,t) = {\{{-ia»{t)a x +1 {t))) , (69) 

and by Wick's theorem one obtains a representation as the Pfaffian of the 2£ x 2£ antisymmetric 
matrix f considered in Ref. [52] 

p?*(e,t) = tf{T), (70) 

which has the same structure of the correlation matrix in Eq. (52) with blocks 

h = ( r ifl -w- 1 ). (7i) 



. ig-i-i ifi 

From the asymptotic correlations (57), we can compute this correlation function in the limit of 
infinite time. Since f%° = 0, the block structure in the matrix T simplifies, and the two-point 
function of the order parameter in the infinite time limit becomes the determinant of a £ x I 
Toeplitz matrix G 

p xx (£,t) = det (G) , Gi n = gfl n . (72) 
The symbol of G is 

3oo(fc) = -e j(fc - efc) cosA fe . (73) 

The representation (72) is more convenient than (70) for determining the large-£ asymptotics of 
p xx (£,t) = det (G). We will show that at late times the two-point function of the order parameter 
decays exponentially with the distance 

p xx {£~> l,t = oo) = C x (£)e~ e ^ [l + o(£ )] . (74) 

The result (74) is obtained by using exact results on the asymptotic behaviour of the determinant 
of Toeplitz matrices, such as Szego's lemma, the Fisher-Hartwig conjecture, and generalizations 
thereof (the results we need in the following are summarized in Appendix A and are taken from 
Refs. [68, 69, 70, 71]). To apply these theorems it is useful to change variable as z — e lk , so that 
the symbol can be rewritten as 

/ \ , x h + h y/z (z- fti)(l//ii - z) 

S 00 (z) = g 00 (-%hiz) = —=r = — , (75) 

2 V' l o y/ho - zy'z - 1/ho z-h 
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where h\ is defined in Eq. (24). The asymptotic behaviour of the determinant and so the physical 
quantities £ and C x (£) depend on the analytic properties of the symbol, which in turn change with 
the quench parameters. 



4-1- Quenches within the ferromagnetic phase. 

The symbol (75) is different from zero on the unit circle \z\ = 1 and has winding number zero about 
the origin z = 0. Under these conditions, it is possible to apply the strong Szego lemma (see (A. 3)). 
The inverse correlation length is straightforwardly obtained as 

r 1 = - ^ln floo (e ifc ) = -j** glncosA fe . (76) 

The integral can be carried out using (x > 0) 
r27T dk 

— \n{x~e lk ) = 6(x-\)\nx, (77) 
to finally obtain 

r^-K^r). (78) 

The multiplicative factor C x (£) = Cf F is slightly more complicated to obtain. It is given by (see 
Eq. (A.3)) 

Cp F = exp^fc(lnfl 00 ) fe (lnfl 00 )_ fe , (79) 

fe>i 

where we used the same notation as in the appendix. The Fourier coefficients (In goo)*; are 

r {h k -2hi k )/(2k) if fc>0, 

(Infloo)fc = I ln((/i + Jio)fci/2) if fc = 0, (80) 

[ {2h\ -2h- k - h- k )/(2k) if k<0. 

Carrying out the fc-sum in (79) we obtain 

C x = l-ft/io + y/(l-fe 2 )(l-fe§) f . 

FF 2VT^hho^/l^hl ' 1 ' 

Putting everything together we conclude that the two-point function has the following asymptotic 
behaviour 

p-(£, t = oo) = 1_ Hh ° + ^^-f^K) e -m + o(e-^) , (82) 

for some k > In Fig. 3 this analytic prediction is compared with numerical data. The 

agreement is clearly excellent. 
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Figure 3. \n\p xx (£,t = oo)| after a quench within the ordered phase from hg = 0.3 to h = 0.7. 
The straight line is the asymptotic prediction in Eq. (82). 



4-2. Quenches within the paramagnetic phase. 

Here the winding number of the symbol equals 1. Following Ref. [71] (see Appendix A for some 
details) we first rewrite the symbol by factorizing an appropriate power of z, such that the remainder 
has winding number zero about the origin 

h + ho \fz (z — l/hi)(hi — z) 



2/0) 



,(l/z)z 



(83) 



2hy/h^ y/ z - 1/hoVha ~z z-l/h 

The function y(z) has winding number zero as desired and no zeroes on the unit circle \z\ = 1. We 
note that ^^(l/z) is in fact the symbol of the transpose of G. It has negative winding number, 
which is precisely the case discussed in Appendix A. From the Wiener-Hopf factorization of y(z) 
as defined in Eq. (A. 9) 



y( z ) = y+( z ) 



h + hr 



hi y-(z) 



2hh 

y±(z) = cxp^2(\n y) ±k z ±k , 

k>l 



the correlation function is obtained using (A. 16) 



; (£, t = oo) ~ exp fc ( ln v)kQn v) 



-k 



k>i 



h + ho V T 



Here Ii is given in terms of y±(z) by 



h = 



2 71 



dfc 

2^ 



i(:k 



; y-(e- ife ) 



o— ik 1 



(84) 



(85) 



(86) 
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The inverse correlation length can be read off from (85) to be 

In order to find y~(z) we use the standard Wiener- Hopf factorizations 

<«<«-*)«)- = { **<}; (88) 

(«(*-*)•)+-{ (w*)« J :>1; ( 8q ) 

where a is an arbitrary complex number. The rhs above are independent on a because of our 
normalization of the factorization in Eq. (A. 9). 

The Wiener-Hopf factorization of y(z) is a combination of the above with q = 0, ±1/2, ±1, and 
after some simple algebra we get 



/ \ z-l/hi y/ho hi-z 

Vz - l/h z ~ l / h h i Vh -z 

Substituting these expressions into It, we obtain 

hVhp r dz e _i y/z-l/fo hl~Z 

H h\ f c 2m Z ^ho—z (h-z)(z-l/hi)' iy±j 

where C is the unit circle. For numerical computations the following equivalent expression is more 
convenient 

h = 2Jh ^-£_ e -*W*)+eo(fc)-2« 1 (fc)) (92) 
J 2ire h (k) 

where 8i is the Bogoliubov angle corresponding to the magnetic field hi, i.e. 

e ie l{ k) = h \ e = . (93) 
\/l + h\ — 2hi cosk 

The integral (91) is dominated by the vicinities of the points in the interior of the unit circle, at 
which the integrand is non-analytic. There are two branch points at z — and z = l/ho, and a 
simple pole at z = 1/hi. If the leading contribution arises from z « l/h±, the ^-dependence takes 
a simple exponential form, while additional power-law corrections arise if the integral is dominated 
by the region z w l/ho- The leading exponential contribution arises from the pole if ho > hi 
(corresponding to ho > h) and by the branch point at z = l/ho if ho < hi (corresponding to 
h > h ) 

T [ h7 l ho > h , 

h a \ , -t > , ( 94 ) 
\ h * h < h . v ' 

Combining this result with (86) and Eq. (87) leads to the following result for the inverse correlation 
length 

r 1 = -ln( / ^/H)+mmin(VM- (95) 
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Figure 4. In \p xx (£, t = oo)| as function of I after two quenches within the disordered phase. Left: 
The case ho > h: The solid line is the prediction (101) and is seen to be in excellent agreement 
with the numerical data. Right: The case ho < h. The straight line is the leading term of the 
prediction (106). Some deviations at small i are visible due to the subleading correction in (106). 



4-2.1. Calculation of the constant Cp P for h > h. From (85) we have 
Cpp = Sexp^fc(lnj/) fe (lny)_ 



)-k 



k>l 



where B is the residue of the integrand in Eq. (91) at l/hi without the exponential factor h 1 
B 



(96) 



(97) 



The Fourier coefficients (lnj/)^ are the same as for the quench from 1/ho to 1/h, because cos Aj, is 
invariant under the inversion of the magnetic fields, so that 



(ln L F ) fc 



hn—th 



(98) 



l/ho-yl/h 



Using this in (96) we conclude that the "Szego part" of Cp P is the same as in the ordered phase 

C FF (l/h ,l/h). (99) 



)-k 



cxp ^/c(lny) fe (lny) 

k>l 

Combining the two pieces in Eq. (96) we have 
C 



•-pp 



lh(ho-h)y/E[=l 

(h + h )(hh -l) ' 



h > h, 



so that the two-point function has the following asymptotic behaviour 

p xx (£,t = oo) 



h(h -h)y/h*-l m U 

(h + h )(hho-l) [ h 



h > h, 



(100) 



(101) 



where k > ^ . The asymptotic result (101) is compared to numerical results for the two-point 
function in Fig. 4 (left panel). The agreement is clearly excellent. 
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4-2.2. Calculation of Cp P (£) for h < h. This case is slightly more involved. The two point 
function acquires a power law prefactor since the leading contribution to the integral (91) arises 
from the vicinity of the branch point at l//io- Isolating this contribution gives 

t -^h- l *- e f 1 ^ ^ ^Z^ZI h.-x/hp un(h -^ (m o\ 

£ ~ /H ° 4 +£l+te 2m ^ ho x/ha (h - x/h )(x/h 1/M + j ' (LU2) 

where e and ei are two infinitesimal positive constants. After the change of variable w = (1 — x)£ s , 
with £ s = £ — | , we obtain 

'~ ^ /I J o 2 ^h -l + I ^Jh-i + ijrM-i-i7t o y 

where 

M, = 4fl-^-ei). (104) 



h 



The large ^ asymptotics is obtained by expanding in powers of l/£ s and taking the integration limit 
(Me) to infinity. At a given order of the asymptotic expansion the resulting error is exponentially 
small. After straightforward algebra we arrive at 

r(3/2) 2hh a Mi-i K l nn ^ 

e ~ 2tt y /f^l(h h-l)(h 1 -h )(£-l)V [ ' 

which in turn leads to the following result for the two-point function 

hoVhfhho - 1 + y/(h 2 -l)(hl-l)) 2 

XXI p . \ ^ V /_ 

9 [ ' >- A^(hi-if/^h h-i)m( h -h ) 

x (£-±y l (l + I ^ T + ...)e-^, h a <h. (106) 

The constant a characterizing the subleading contribution in (106) is given by 
3 4 h (5h (h - h ) - 8^(h 2 - l)(hl - 1))- 



Vhl-1 h h-l (h - h Q )(h 2 Q - 1) (107) 

The requirement that the term involving a should be small establishes the regime of validity of the 
expansion (106). In particular, for small quenches in the sense of paper I, i.e. r\ = max | sin Afc| 
being small, we have 

Hence (106) holds only for distances that are sufficiently large compared to (h—ho)^ 1 . In particular, 
in the limit h — > ho its region of validity disappears entirely. 

In Fig. 4 the asymptotic result (106) (with a set to zero) is compared to numerical data. It is 
evident that for h < h and £ not large enough there are sizeable corrections to (106). 
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Figure 5. In \ p xx (£, t = oo)| as a function of I for quenches across the critical point. Blue squares 
correspond to a quench from the paramagnetic (ho = 1.3) to the ferromagnetic (h = 0.7) phase. 
The blue solid line is the asymptotic expression (118). Red triangles correspond to a quench from 
the ferromagnetic (ho = 0.7) to the paramagnetic (h = 1.3) phase. The red straight line is the 
asymptotic result (116). In both cases the asymptotic and numerical results are seen to be in 
excellent agreement. 



4-3. Quenches across the critical point. 

For quenches across the critical point hi is a complex number with unit modulus \hi\ = 1. In this 
case, the symbol has two simple zeroes on the integration contour at z\ — h\ and z% = 1/hi = h*. 
Thus, as first step we need to cast the symbol in the form (A. 7) which is suitable for application of 
the Fisher-Hartwig conjecture as reported in (A. 10). Since the symbol has no discontinuities, we 
only need to isolate the zeroes. This can be done by rewriting goo(z) in the form 



0? 



(z) = Ih-^ihl)-^} fM\z) z ft+& (-l)A-A (fel ~ z)(z 



K) 



where 



/ 



(/3i,/3 2 ; 



h^(hlf- 



h + h 



(-1)* 



/3i-/3 2 



Vho~ 



The symbol (109) is now in the desired form (A. 7) with 

R = 2 , a = f3 = , a x = - , fa = - - m x , 



1 

0-1 = - 

2 



02 



1 

2 -m 2 . 



(109) 



(110) 



(111) 



where mi. 2 € Z are integers that label the inequivalent representations of the symbol. Taking into 
account the contributions of all representations results in an expression for the two-point function 
of the form 



p xx (l > l,t 



00 



E 

{/3i,/3 2 }i, 



*\02 



-t/C(Pl,fc) 



(112) 
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where we isolated the oscillatory factor [ — h^ 1 (hl)^ 2 ^Y . In (112), £(/?i,/?2) and Cp lt p 2 are 
respectively the correlation length and correlation amplitude corresponding to the representation 
labelled by of the symbol (c/. (A. 10)). In the special case R = 2, ot\ = a 2 = 1/2, a = j3 = 0, 
z\ = hi, Z2 = h\, the amplitudes of the various representations are given by 



9>1 



fi^\h,) 



(A A) 



x [if lA) (M) 
1 11 Khxe™) 



?2_il 2 



&)G(§ 



ft) 



G(2) 



(113) 



where G is Barnes' G-function and f^ 1 '^ 2 \z) are the factors of the Wiener-Hopf factorization of 
j(hJh)^ m ^ e conventions (A. 9). As G(—n) = for n = 0, 1, 2 . . . we need to consider only the 
four inequivalent representations 

a =4- 132 =± \- (n4) 

In order to obtain explicit expressions for C±i ± i and £(±|,±|) we need to consider quenches 
from the paramagnetic to the ferromagnetic phase and vice versa separately. This is done in the 
following two subsections. 



4-3.1. Quenches from the paramagnetic to the ferromagnetic phase. The leading contributions in 
the large-^ limit arise from the two representations 



In both cases the inverse correlation length takes the simple form 

/■27T 



1 1 

±2' T 2 



< ^(±- -) = -jr'gln(^f/(±i.=Fi)(e tt ))=-ln( 
As shown below, the constants Cp lt p 2 are equal to 

1 I ho — h 
2 



2ft 



Gi 



G_i i 

2 ' 2 



(116) 



(117) 



Putting everything together we arrive at the following result for the asymptotics of the two-point 



function 



P xx (e,t 



oo 



ho — h 



1 



Re[h{] e- £ /« + Q{e- K ' 1 ) , 



(118) 



where k? > Rewriting hi as a function of h and ho one recovers the form given in (27). The 
expression (118) is compared to numerical results for p xx (£,t — oo) in Fig. 5. The agreement for 
large distances £ is clearly excellent. We note that (118) can be understood as the leading term in 
an asymptotic expansion only if Re[/if] is of order 1. For the particular quench from 

h + V2{l-h 2 ) 



ho 



(119) 



1 - 2ft 2 

to ft, < l/y/2 we have h\ = i and (118) vanishes for I = 4n + 2, with n integer. For this particular 
sequence of distances £ does then not represent the correlation length. 
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Calculation of the constant Ci _i. The Wiener-Hopf factorization of f^~^(z) is 

/^■-»)W = /i i '- i) W^V»i^'-' ) W ) (120) 

where 

Vho-z y/z-l/h a z-h 
This results in Fourier coefficients 

( ho k /{2k) if £;> , 

(ln/(5.-D) fc = J In(^fti) if fc = , (122) 

[ -{h* + 2h- k )/(2k) if fc<0. 

The Szego part of the constant is then of the form 

-[E/(W^-\(W^-»)_J = ^=. (123) 

Substituting (121) and (123) into the expression (113) we obtain the result (117). The calculation 
of C_i i is completely analogous. 

4-3.2. Quenches from the ordered to the disordered phase. Here the leading contribution in (112) 
arises from the representation 

0i=ft = ^. (124) 
The Wiener-Hopf factorization of f^'^(z) is 



where 



f^\z) = f^(z)^f^\z), (125) 



tf'» ) (z) = 7f i Tr ^-, f M\ z)= ^*. (126) 
VI — z/i 1 — z/h v'z — ho 

The Fourier coefficients of In fd' i~> are 

r (ftg +2/i- fe )/(2/c) if k > , 

(ln.f(^)) fe= l n (M^) if fc = 0, (127) 

[ -h^ k /(2k) if fc<0. 
while the inverse correlation length is 

r 1 = (ln/^) =ln(^). (128) 

We note that (128) coincides with the result obtained by formally interchanging ho with h for 
quenches from the disordered to the ordered phases (116). The Szego part of the constant is given 

by 

-[EMM^^(M'»lJ= (/p | p =. (129) 



k>l 
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Substituting (126) and (129) into (113) and (126) we obtain the following result for the asymptotic 
behaviour of the two-point function 



where k > £ 1 . In Fig. 5 the analytic result (130) is compared to a numerical computation of 
p xx (£, t = oo) and the agreement is seen to be excellent. 

4-4- A general formula for the correlation length. 

The results (78), (95), (116) and (128) can be combined into the following general expression for 
the inverse correlation length 



where h\ is defined in Eq. (24). After simple algebraic manipulations, this expression can be cast 
in the form reported in (23). 

5. Transverse correlations 

We now turn to the correlation functions of transverse spins of. These are simple because of are 
local in the fermionic representation of = ia% af . Concomitantly, and in contrast to order parameter 
correlators, exact analytic results have been known for many decades, both in equilibrium [72, 52] 
and after a quantum quench [52]. Here we summarize these results and point out some interesting 
features, which as far as we know have not previously been stressed in the literature. 

5.1. The one-point function 

We first consider the one-point function, i.e. the magnetization in the z direction, which is [52] 

f n dk / \ 
(erf) = (iaf {t)af (t)) = - J — [cos 6 k cos A k + sm9 k sin A k cos(2e h (k)t)J , (132) 

where Sh{k), 9 k , A k and defined in (3), (10) and (11) respectively. Unlike the order parameter, the 
transverse magnetization has a finite value in the stationary state 



The result (133) agrees with the GGE prediction. The late-time limit of (132) can be determined 
by a stationary phase approximation. The leading contributions arise from the saddle points at 
k = and k — tt and a straightforward calculation gives 




(130) 




(131) 




(133) 




(134) 



+ 



h -h rsin(4J|l - h\t + tt/4) sin(4J(l + h)t - tt/4) 



+ 0((Jt)-S). (135) 



4^{2hJt)i L v/fl^lll-Ziol VT+h(l + h ) 
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5.2. The two-point function 

The connected transverse two-point correlator is expressed in terms of the Majorana fermions as 
(seeEq. (53)) 

= «4fi) « a ?+i> - K a ?+i) « a £+i> = f-ife ~ 9e+i9i 
Using Eq. (55), p zz (£,t) can be written as a double integral in momentum space 

2 



pz ^ t] = £ w e " (fel " fe2) { n sinA ^ in (2^)*) 



3 
2 

n 



i6 k . 

e •» 



cos A kj - isin A fcj cos (2e h (kj)t) 



where 8k, A^ and Eh(k) are given by (10), (11) and (3) respectively. 



(136) 



(137) 



5.2.1. The infinite time limit and the GGE. In the limit t — »■ oo at fixed, finite t all oscillating 
terms in this integral average to zero and we are left with 



pT(*,oo) = 



^ e uk e ie k 
2ir 



cos A fc 



dp. 



e «P e -ie p cosA 



2tT - p 

This resembles the result for the finite temperature connected 2-point function obtained in [52 

dfc 



(138) 



z\\2 



«"fof + /» " «*f » 



/* ge^e^tanh(^M) |* ge^e"*' tanh(^) . (139) 

Here ((C)) denotes the thermal equilibrium expectation value at temperature l/f3. Replacing j3 
with a mode dependent inverse temperature as in (65), 

(3e h (k) — > (3 k e h (k) = 2 arctanh( cos(A fe )) , (140) 

reduces (139) to (138). This again confirms the notion that the GGE can be thought of in terms 
of a mode-dependent temperature. In order to determine the large— I behaviour of p zz (£,oo) we 
change the integration variable in (138) to z = e lk . This results in a contour integral along the unit 
circle C 

(h + h f 
4hh + ' 
dz z £ - 1±1/2 (z-l/h^h.-z) 

2™ y/( z - l/h )(h - z) ±(z-h* 1 ) ' 

where h\ is given in (24). Inside the unit circle there is a branch cut connecting with min{/io, 1 /ho}. 
In addition there is a pole at mm{h, 1 /h} in I + if h > 1 and in J_ if h < 1. The leading contributions 
to the integral arise from the pole and the vicinity of the branch point at min{ft,o> 1/ho}. Taking 
these into account (and assuming that hgh ^ 1) we obtain 



/Woo) 



-2\\nh \l 

cr 5 — = — (i 



or 1 )) 



if | lnh\ > | \nh \ , 



-(|lnh| + |lnh |)€ 
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-(l + Oir 1 )) if |In/io| > |ln/i|. 



(142) 
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Here the amplitudes are 

\h -l/h \ h 



C{ z = 



Air 1 - hh 



(h-l/h)y/\ho-l/ho\(ho-h) 



oS gn(ln/i)|ln h \/2 



2 8y^h \] h (hh -l) s - mh \^h\ + \lnho\ ■ y °) 

For quenches such that hho = 1 the result is instead 

pf(£,oo) = - (/l ~ 1/fe)2 e- 2 l 1 " fe l £ + ... h = l/h. (144) 

Z7T 

5.2.2. Time dependence in the space-time scaling regime. We now consider the late time behaviour 
in the space-time scaling regime, i.e. in an asymptotic expansion around the limit £, t — > oo at 
fixed finite ratio v = £/2t. As I oc t appears in (137) linearly in the phases (albeit in different 
combinations), the asymptotic power law behaviour can be determined by a stationary phase 
approximation. The calculation is relatively simple, but the resulting formulas are long and not 
very illuminating. The general form of the answer is 

p*'(£, t)-]j2 a ^ v ) cosM«)t + <pi(v)) , (145) 

i 

where the sum runs over the various saddle points. The leading power- law i _1 comes simply 
from the stationary phase, but the calculations of the functions a,, ipi and u>i are quite tedious. 
For this reason we focus on the non-oscillating contribution with uiq(v) = 0, which we denote by 
PQ Z (£,t). This zero-frequency piece arises from terms, in which the phase in Eq. (137) vanishes at 
the stationary points, i.e. from 



/ 



d|^2 e 2i„t( fcl - fca ) sin Afci sin Afc2 cos (2Mfci) _ £h{k2)]£) (1 + e i(e kl+ e k2)) (M6) 



The large-t behaviour of this double integral can be extracted by a two-dimensional stationary 
phase approximation. Its non-oscillating part (characterized by k\ — k 2 ) for h ^ 1 is 

pl*{£, t)c±-l — sin 2 A k cos 2 k S(e' h (k) - v) (147) 
4i Jo 7r 



?g("max-«) r sin 2 A fcv cos 2 K sin 2 cos 2 9^ 



W'(k v )\ \e"(k v )' (J4<S) 

where k v and k v are the solutions of s' h (k) = v. In the case h = 1 there is only a single term in 
(148) with k v — 2arccos(u/2 J). The following remarks are in order 

(i) For v > u max the connected two-point function is exponentially small in vt = £ » 1, in 
accordance with expectations based on causality [11]. 

(ii) In contrast to the order parameter two-point function p* x (£ 7 t) studied in paper I, for given t 
and £ only elementary excitations with velocity exactly equal to v = ±£/2t contribute to the 
asymptotic behaviour of p z c z {£, t) . This is a consequence of the transverse spins being fermion 
bilinears. 

(iii) In the limit v — > (147) becomes (sin A k oc v, while cos 2 9 k and e'^(k) are finite at the saddle 
points), 

P z z (£,t) <x£ 2 /t 3 . (149) 
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6. Conclusions 

In this work we have derived analytic expressions for the stationary behaviour of the reduced density 
matrix as well as one and two point functions in the transverse field Ising chain after a sudden quench 
of the magnetic field. We have shown that local observables in the infinite time limit are described by 
an appropriately defined generalized Gibbs ensemble. We have furthermore analyzed the approach 
in time to the stationary state and addressed the question, how long it takes for the GGE to reveal 
itself for a given observable. We found that in general this time grows exponentially with the spatial 
extent of the observable considered, e.g. t oc e ae for a two-point function with spatial separation £. 
Our analysis can be straightforwardly generalized to quantum quenches in other models with free 
fermionic spectrum such as the spin-1/2 XY chain in a magnetic field and we expect qualitatively 
similar behaviour. In contrast, the stationary behaviour in interacting integrable models such as 
the (5-function Bose gas is not amenable to treatment by the methods employed here and one needs 
to proceed along different paths [42]. 
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Appendix A. Asymptotic behaviour of determinants of Toeplitz matrix: Szego 
Lemma, Fisher Hartwig conjecture and generalization 

In this appendix we summarize results on Toeplitz determinants, that are used in the main part of 
the paper. Let (Tg) = ti_ n be a I x I Toeplitz matrix with symbol t{e lk ), i.e. 



The behaviour of the determinant of T for asymptotically large I depends on the analytic properties 
of the symbol t(e tk ). 

Appendix A. 1. Smooth symbol: The strong Szego Lemma 

If the symbol t(z) is a nonzero continuous function of z on the integration countour, with winding 
number zero about the origin, the function logi(z) admits the Fourier (Laurent) expansion 



The asymptotic behaviour of the determinant is then given by the strong Szego limit theorem [68] , 
which states that 




(A.l) 





(A.2) 



det [T/] = E [t] e ^ lost )°(l + Oil 1 - 213 )) 



(A.3) 
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(A.4) 
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where 

E [t ] =exp[^g(logt),(logi)_ 

A convenient alternative form of this result is 

In det [2i] / 27r ^lnt( e ifc ) + V g (logt) 9 (logt)_ 9 + 0(^-^). (A.5) 

The exponent (3 characterizing the subleading corrections is determined by the analytic properties 
of t(e lk ). The integer part of (3 is equal to the number of continuous derivatives of t(e lk ). If the 
symbol is infinitely differentiable the corrections fall off faster than any power of I. 

Appendix A. 2. Symbols with zeroes and singularities but zero winding number: The Fisher 
Hartwig conjecture 

One of the many generalizations of the Szego theorem is the Fisher Hartwig conjecture. It applies 
to cases in which the symbol t(e ik ) has zeroes and/or discontinuities and can be expressed in the 
form 

R 

t(e lk ) = t a {e lk ) J[ e *A-(*-fcr-iragn(k-fer))(2 - 2cos(fc - k r )) a " . (A.6) 

T-0 

Here R is an integer, a ri f3 ri and = k < ki < ■ ■ ■ < < 2ir are constants characterizing the 
location and nature of singularities and zeroes, and to(e lk ) is a smooth nonvanishing function with 
winding number zero. In terms of the variable z = e tk the symbol is 

*(*) = toWzVU* fl( 2 - 7 - fJW*)^ ■ ( AJ ) 

3=0 3 

where Zj = e lkj , j = 0, . . . , R and 

e -A- 0<argz<fc J 

e~™^ kj < arg z < 2tt. [ ' 

We follow conventions in which z = 1 is always included in the set {zj}, while for all other Zj we 
must have either ctj ^ 0, f3j ^ or both. Following the notations of Ref. [69] we introduce the 
Wiencr-Hopf factorization of t (z) as 

t (z) = 6 + (z)e( losto )°6_(z), 

b+(z) = e ^->^ ost ^ z \ b-(z) = e ^^-i( lo s*)^ 9 . (A.9) 

Here the constant piece e' 108 * -* has been fixed through the requirement (lnfe + )o = 0. When 
Re[aj] > —1/2 and for aj ± f3j £ 7L~ , the Fisher-Hartwig conjecture (which under the above 
conditions has been proved [70]) gives the asymptotic behaviour of the determinant 

l[to] e *(io gto )o M*i)r°" +ft M*)]-"*) FUtf-ft 



9p 3 ( z ) 



det 



E\ 



n 



0<j<k<R. 
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n 

3=0 



G(l 



ft) 



(A.10) 



Here the functional Em is given in (A. 3) and G(x) is the Barnes G function. In the special case 
R = 0, a = /3 = the strong Szego lemma applies and the above formula reduces to (A. 3). In 
cases where ctj and 0j don't satisfy the above requirements, the large-£ asymptotic behaviour can 
be determined from the so-called generalized Fisher-Hartwig conjecture. Now there are in general 
several different representations of the symbol in the form (A. 6), i.e. there are several possible 
choices of the parameters {(3 r } in (A. 6). The asymptotic behaviour of the determinant is then 
given as a sum over all inequivalent representations 

det[2>]~ ]T detpi] {/Jp}) (A.ll) 

incq 

where det \Tf\ip j denotes the expression (A. 10) for a given set {f3 r }. In the cases encountered in 
the main part of our paper we identify the representations giving rise to the leading asymptotics 
and quote only their contribution. 

Appendix A. 3. The case of a symbol with non-zero winding number 

A generalization of the Szego theorem for symbols with nonzero winding number does exist [71] . If 
the symbol has negative winding number, i. e. 3k e N + such that 

a(e lk ) = (-l) K e lKk t{e lk ) (A.12) 

has winding number zero about the origin, the large-£ asymptotics of the determinant is given by 
[68, 69, 71] 



dot [T e ] = E [a] exp(tj 
V.3) a 



dk 

2^ 



\n[a{e lk )} det 



0(r 3/3 ))(l + 0(^- 2/3 )) , (A.13) 



where E^ is given by (A. 3) and T K is a k x k Toeplitz matrix with elements 

* r -i(l-n)k -Uk a -( e% ) 



2tt 



a + (e ik ) 



where 



oo 

a(e ik ) = a + {e lk )e^ osa ^a-{e lk ) a ± (e lk ) = cxp[^(lna) ±J 



,±ijk 



(A.14) 



(A.15) 



As before we use notations where (lna)j are the Fourier coefficients of In a. For our purposes we 
only need to consider the case k = 1, in which (A.13) takes the simpler form 



detT e ~E [a] e 



dfc 
27' 



a_{e ik ) 
a + (e lk ) 



(A.16) 



dfc 



-ilk 



a_(e 4fe ) 



(A.17) 



At large I we have det[T^] oc e e ^, where 

r-i i- detT e n ^ 1, f 2 

t = hm — - — - = (loga)n + urn - In / 

The case of a symbol with positive winding number can be obtained directly from (A.13) and 
(A.16) by noting that the transpose of a Toeplitz matrix is another Toeplitz matrix with a symbol 
of opposite winding number. 
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